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Abstract —The AC power flow equations underlie all opera¬ 
tional aspects of power systems. They are solved routinely in 
operational practice using the Newton-Raphson method and its 
variants. These methods work well given a good initial “guess” 
for the solution, which is always available in normal system 
operations. However, with the increase in levels of intermittent 
generation, the assumption of a good initial guess always being 
available is no longer valid. In this paper, we solve this problem 
using the theory of monotone operators. We show that it is 
possible to compute (using an offline optimization) a “mono¬ 
tonicity domain” in the space of voltage phasors. Given this 
domain, there is a simple efficient algorithm that will either And 
a solution in the domain, or provably certify that no solutions 
exist in it. We validate the approach on several IEEE test cases 
and demonstrate that the offline optimization can be performed 
tractably and the computed “monotonicity domain” includes all 
practically relevant power flow solutions. 


1. Introduction 

Power systems are experiencing revolutionary changes due 
to various factors, including integration of renewable gener¬ 
ation, distributed generation, smart metering, direct or price- 
based load-control in operations. While potentially contribut¬ 
ing to the long-term sustainability of the power grid, these 
developments also pose significant operational challenges by 
making the power system inherently stochastic and inhomoge¬ 
neous. As these changes become more widespread, the system 
operators will no longer have the luxury of large positive and 
negative reserves. Thus, operating the future power grid will 
require developing new computational tools that can assess the 
system state and security margins more accurately and faster 
than current approaches. Specifically, these new techniques 
need to go beyond linear models and ensure that the power 
system is safe even in the presence of large disturbances and 
uncertainty, where nonlinear effects dominate. In this paper, 
we focus on the fundamental equations of the power system : 
the power flow (PF) equations. The PF equations constitute 
a system of nonlinear equations and are known to exhibit 
complex and chaotic behavior 0 0 - Standard techniques like 
Newton-Raphson and its variants often fail to converge when 
the operating conditions are changing rapidly or the system is 
close to its security margins. In such a situation, it becomes 
difficult to assess whether the system is actually operationally 
unsafe or if the Newton-Raphson method failed because of 
numerical difficulties or bad initialization. In this paper, we 
propose an approach to remedy this problem. Our approach 


is based on the theory of monotone operators Just as a 
convex optimization problem can be solved efficiently, one 
can find zeros of a monotone operator efficiently (in fact, 
the former is a particular case of the latter as the gradient 
of a convex function is a monotone operator.) Thus, if we can 
show that the nonlinear PF equations can be described by a 
monotone operator over a sufficiently large domain, then they 
can be solved within the domain of monotonicity efficiently, 
or certify that no solution exists within the domain. . It turns 
out that the PF operator is not globally monotone, however it 
is monotone over a restricted domain. 

Our main contribution is an efficient procedure based on 
semidefinite programming to characterize an operationally 
relevant domain (that contains all the voltages satisfying 
operational constraints) in the space of voltages over which 
the power fiow operator is monotone. The domain is specified 
in terms of a simple bound on the voltage phasor ratios 
between neighboring buses. Roughly speaking, these bounds 
require that the voltage phasors between neighboring buses are 
“sufficiently close”. Once this is done, there are well-known 
efficient algorithms to compute solutions of the PF equations 
satisfying these bounds, or certify that no solutions exist 
within the monotonicity domain. These algorithms are based 
on solving an associated monotone variational inequality, for 
which several theoretically and practically efficient algorithms 
have been developed (see and chapter 12). 

As a by-product of our analysis, we obtain a domain over 
which the power fiow Jacobian is non-singular. The domain is 
characterized by simple bounds on the voltage phasor ratios 
between neighboring buses. The bounds are interesting in 
their own right, as they allow the system operator to monitor 
distance to voltage collapse or loss of synchrony (where 
the Jacobian becomes singular) simply by monitoring ratios 
between voltage phasors on neighboring lines, which can even 
be done in a distributed manner. Numerical tests show that the 
bounds obtained are non-conservative and cover a wide range 
of operating conditions. 

The rest of this paper is organized as follows. Section 
covers relevant background on power systems and monotone 
operators. The main technical results are presented in Section 


|IIl| In Section |IVj we discuss how our approach compares 
to related work. In Section [Vj we present numerical results 
illustrating our approach on some IEEE benchmark networks. 


II. Modeling Power Systems 


K.D. and S.L are with the department of Computing and Mathematical 
Sciences, Caltech, M.C is with the T-Division, Los Alamos National Labora¬ 
tory, Los Alamos, NM, 87544 USA. All queries should be addressed to K.D 
(e-mail: dvij@caltech.edu). 


A. Notation 

M is the set of real numbers, C the set of complex num¬ 
bers. denote the corresponding Euclidean space in n 
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dimensions. Given a set P C Int {V) denotes the interior 
of the set. Given a complex number x G C, Re {x) denotes 
its real part and Im (x) its imaginary part, x* denotes its 
complex conjugate. \\x\\ refers to the Euclidean norm of a 
vector X G or X G and (x, y) to the standard Euclidean 
dot product. Given an vector x G di {x) denotes the n x n 
diagonal matrix with (i,i)-th entry equal to Xi.The nuclear 
norm of M is denoted by and is equal to the sum of its 

singular values. Given a differentiable function / : i-G 

V/ denotes the Jacobian of f, Sikxk matrix with the i-th row 
being the gradient of the i-th component of /. Eor M G 
Sy (M) = Given two sets of indices S, S' and matrix 

M whose rows and columns are indexed by 5", = [M]^ 

denotes the |5'| x |5'| matrix with the following property: 

Mij ifijeSnS' 

0 otherwise 


exp(p+j6>). We also write V = exp^ and whenever 

we write V it is understood that the constraints on | Vp^ \, Vq 
are satisfied. Let Oij = Oi — Oj , pij = Pi — Pj. 

Definition 1 (Power Elow Operator). Define the power fiow 
operator F as 

n 

[F {V)]i = Bij exp {pi + Pj) sin (%) 

j=0 

n 

+ F “ PiT < * < «• (la) 

i=o 

n 

[F{V)Ui = Y. Gij exp {pi + Pj) sin (%) 

j=0 

n 

- yy Bij exp {pi + Pj) cos {Oij) -qi,l<i< |pq| (lb) 

i=o 


tril (M) denotes the vector formed by the lower triangular 
entries of matrix M. 

B. Background 

We represent the transmission network as a graph (V, f) 
where V is the set of nodes and £ is the set of edges. In 
power systems terminology, the nodes represent the buses and 
the edges correpond to power lines. Buses are denoted by 
indices i = 0,1,... ,n and lines by ordered pairs of nodes 
We pick an arbitrary orientation for each edge, so that 
for an edge between i and j, only one of (i,j) and (j, i) is 
in £. If there is an edge between buses i and j, we write 
i jj ^ i. 

The transmission network is characterized by its complex 
admittance matrix Y G F is symmetric but not 

necessarily Hermitian. Let G = Re (F), B = lm (Y). 

Let Vi be the voltage phasor, pi and qi denote active 
and reactive injection at the bus i respectively. V is the 
vector of voltage phasors at all buses. Three types of buses 
are considered in this work: PV buses where active power 
injection and voltage magnitude are fixed, while voltage phase 
and reactive power are variables. The set of PV buses is 
denoted by pv. The voltage magnitude set point at bus i G pv 
is denoted by Vi. 

PQ buses where active and reactive power injections are fixed, 
while voltage phase and magnitude are variables. The set of 
PQ buses is denoted by pq. 

Slack bus , a reference bus at which the voltage magnitude and 
phase are fixed, and the active and reactive power injections 
are free variables. We choose bus 0 as the slack bus as a 
convention. The slack bus is voltage phasor by Vq. 

We denote the union of PV and PQ buses as nsb = pv U pq. 

We will work with the logarithmic-polar representation of 
the voltage phasor: 

Vi = exp {pi + j0i) , Pi = log (IViI), 6»i = ZV^i 

The variables in the power fiow problem are the phases at 
the non-slack buses O^sh and the voltage magnitudes at the 

PQ buses ppq. Eor brevity we write where V = 


Remark 1. This is a simplified version of the practical power 
fiow problem: We do not account for transformers (considering 
the entire network in renormalized voltage units) and we 
also do not consider more general 7r-model for power lines, 
accounting for shunt capacitors to the ground. However, all the 
aforementioned (and other) important practical details can be 
easily incorporated in the model we use and are not restrictive 
in terms of applications of our results to practical power 
systems. 

We denote by pqq = pq + n the set of indices of the 
PQ buses shifted by n. Thus, the power fiow operator F 
is indexed by nsb U pqq. The power fiow equations can 
be written as F {V) = 0 solved for denote by 

Jf {V) the Jacobian of F with respect to evaluated at 
V = exp^ Note that this is not the standard power 

fiow Jacobian in the power systems, since we differentiate wrt 
log (I V|pq) rather than | V|pq. We denote by k the total number 
of variables being solved for (k = |nsb| + |pq|). 

The next lemma expresses Jf {V) as a quadratic function 
of V: 


lemma 1. The power flow Jacobian Jf {V) G can he 

written as a quadratic matrix function of the voltage phasors: 

Y,^i\Vp+ Y ^ijBe{ViVj*) + <i>ijlm{ViVj*) (2) 


where 


A = 


= 


V — 

V ij — 


\f 
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Proof: Via direct differentiation. ■ 

Remark 2. The Jacobian need not be symmetric, as can be seen 
by the lack of symmetry in the matrices A^, However, 

each entry of the Jacobian matrix is a quadratic function of 
the voltage phasor V. 

C. Monotone Operators 

We now review briefly the theory of monotone operators, as 
is relevant to the approach developed in this paper. For details 
and proofs of the results quoted in this section, we refer the 
reader to the recent survey j^. A function Ff : i-A is 

said to be a monotone operator over a convex domain V if 

{H (x) - H {y),x - y) >0 'ix,y eV 

A monotone operator is a generalization of a monotonically 
increasing function (indeed, if k = 1, the above condition 
is equivalent to monotone increase: x > y H (x) > 

H (y)). H is said to be strictly monotone over V if 

{H {x) - H {y),x-y) > 0 \/x,y e T>,x y 

A common example of a monotone operator is the gradient of 
a differentiable convex function. 

Definition 2 (Monotone Variational Inequality). Let P c 
be a convex set and H he sl monotone operator over V. The 
variational inequality (VI) problem associated with H and V 
is: 

Find X eV such that {H {x),y — x) >0 My e V (6) 

Deflne the normal cone to P at x as Mv {x) = {y : 
(^, z — x) < QNz G V}. Then, the variational inequality is 
equivalent to flnding x e V such that —H{x) G A/i) (x). 

The following result shows that monotone variational in¬ 
equalities with compact domains always have a solution and 
can be solved efficiently. 



(a) VI: solution in interior 



III. Monotonicity of the Power Flow Operator 

In this section, we study the monotonicity of the PF operator 
As described in Section |II-C[ zeros of F (solutions to 
the PF equations) can be found efficiently if F is monotone. 
Thus, if we can prove that the PF operator is monotone, the 
PF solutions can be found efficiently. Since PF equations can 
have multiple isolated solutions, it is not possible that the 
PF operator is globally monotone because this would imply 
a unique solution to the PF equations. Thus, we focus on 
characterizing domains over which the PF operator (or a scaled 
version of it) is monotone. This leads to the constrained power 
flow problem: 


Theorem II.l. If H is strictly monotone operator over a 
compact domain V, then (|^ has a unique solution x*. Further, 
an approximate solution x^ G V satisfying 

||a:;e - a:|| < e (7) 

can be found using at most O (log (i)) evaluations of H and 
projections onto V. 


Remark 3. In this manuscript, we are interested in finding 
zeros of the PF operators introduced above. We can use 
monotone operator theory for this as follows: Suppose H 
satisfies the hypotheses of theorem |II.1| If there exists a point 
X* G P with H (x*) = 0, then this is the unique solution of the 
variational inequality (figure Conversely, if the variational 
inequality has a solution with H (x*) 7 ^ 0 (figure lb), then 
have a certificate that there is no solution of Ff (x) = 0 with 
X eV. 

The next result provides a simple characterization of mono¬ 
tonicity for differentiable operators: 


Theorem II.2. Suppose that H is differentiable. Then H is 
strictly monotone over V if and only if 

WH{x) + WH{xf yO WxgV 


Definition 3 (Constrained Power Flow Problem). Given a set 
P C the constrained power flow problem is to determine 
whether the power flow equations F {V) = 0 have a solution 
with G V, and if so, compute the solution. We denote 
this problem by PF (V). 

Remark 4. In this paper, we will characterize sets V such that 
the constrained power flow problem can be solved. 


A. Characterization of Domains of Monotonicity of the Power 
Flow Operator 

We now derive a procedure to characterize the domain of 
monotonicity of the PF operator Q- We first note that solving 
the equations F (V) = 0 is equivalent to solving the equations 
WF{V) = 0 for an invertible matrix W G Define 

Fw (x) = WF (exp^ (x)) for brevity. 

Theorem III.l. Let V be any compact convex set in such 
that 3W G satisfying 

Sy {WJf {V)) yO MV : G V ( 8 ) 

Then, Fw is strictly monotone over the set V. Let x"^ be the 
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unique solution to the monotone variational inequality: 

Find X e V 

{Fw {x),y-x)>0 \/y eV (9a) 

Then, if F (exp'^ (a;*)) = 0, exp'^ (a;*) is the unique solution 
of the PF equations in V. Otherwise, there are no solutions 
of F (y) = 0, G V. Thus, the constrained power flow 
problem can solved efficiently with domain V. 

Proof: This is a straightforward application of the theory 
of monotone operators. A similar proof can be found in 0 - 

■ 

Remark 5. Theorem m characterizes condition 0 under 
which the constrained power flow problem can be solved. In 
the reminder of this section, we show how one can construct 
sets V algorithmically to satisfy this condition. Notice that 0 
is an instance of a containment problem that has been studied 
in the optimization literature j^. 

B. Monotonicity Domains: 2-bus Network 

In order to motivate our choice of V, we first consider a 
simple 2-bus network. Bus 0 is the slack bus and bus 1 is 
a PQ bus. Let y = g — denote the complex admittance 
(conductance g, susceptance h) of the line between 0 and 1 
and ignore any shunt elements. Let the slack bus voltage be 
Vb = 1 (magnitude 1 , zero phase) and the voltage phasor at 
bus 1 be exp {p -y]0). The PF equations are given by 

Pi = 6 exp (p) sin {0) — g exp (p) cos {0) + g exp ( 2 p) 
qi = —pexp (p) sin {0) — 6 exp (p) cos {0) + 6 exp ( 2 p) 

The power flow Jacobian (scaled by exp (—p)) is given by 

b cos {0) + g sin (0) 2g exp (p) — g cos {0) + b sin {0)\ 
—g cos (o) + b sin (<9) 2b exp (p) — g sin (0) — b cos (0) J 

Choosing W = 
comes 

. . . . f cos {0) sin ( 6 >) A 

^ ^ ^ ^ ysm (6) 2 exp (p) — cos [6) J 

The condition ^yiWJp (p,^)) ^ 0 reduces to the diagonal 
entries and the determinant being positive, which simplifies 
to: 

cos ( 6 ») > i exp (-p) = \^y ( 10 ) 

This is a well-known voltage stability criterion for the two- 
bus network Q. In this case, it is also easy to see that when 
cos (P) = the Jacobian is in fact singular, so that this 

is “maximal ’ monotonicity domain, in the sense that there is 
no monotonicity domain that contains it. In fact, the “high 
voltage” branch of the power flow solution set will lie always 
within this domain, until the point of maximum loadability is 
reached (beyond this point there are no solutions to the PF 
equations). 

This example also shows that the choice of W is critical. 
For example, if we choose W = I, then we would automat¬ 
ically require b cos (0) g sin (0) > 0 (in addition to other 
conditions),leading to more restrictive constraints than ( p^ . 


C. Algorithmic Computation of the Monotonicity Domain 

For general networks, the situation is not as simple as in the 
case of the 2-bus network. There is no simple analytical choice 
of W that determines a large monotonicity domain. Instead, 
we propose a computational technique based on semideflnite 
programming. Motivated by the form of the 2-bus constraint 
we consider domain V (7) defined as: 

{(('nsb,/3pq) : COS (ffij) > exp {\pij\) e £} (11a) 

= {: Re max (1^7, ) V (i, j) e £:} 

(11b) 

In the special case of the 2-bus network, the condition ( pTj ) 
reduces to cos{0) > 7exp(p), and as we saw before, one 
can choose 7 = ^ • The specific form of (7) is convenient 
for our purposes and produces non-conservative estimates of 
the true monotonicity domain. Further, the condition reduces 
to simple bounds on the voltage phasors at the end of each 
transmission line, and may be useful in stability monitoring. 
However, depending on particular operational constraints at 
play in the system, we may consider other domains as well. 

D. Computation of the Monotonicity Domain 

Certifying that V (7) is a monotonicity domain (i.e, check¬ 
ing that it satisfies condition ( p^ ) amounts to checking that 
the following problem is feasible: 


Find W such that 

Sy {WJf {V)) ^0 VR ; e V ( 12 a) 

Rewriting this explicitly, we need to solve the following 
optimization problem: 

max min z^Sy (WJp (V)) z (13a) 

w zeR^yec^ 

Subject to Re (R^R/) > yij max {\vy, ) (13b) 

\vy = G pv (13c) 

Ro = 1 (13d) 

z^z = 1 (13e) 

lll^llnuc < 1 (13f) 


The constraint ||lL||nuc — 1 is an arbitrary scaling constraint 
(since the problem is homogeneous in W) where H-H^ denotes 
the nuclear norm, or the sum of singular values of W. If the 
optimal value of the above problem is positive, the optimal 
solution IL* is such that Fw* is strictly monotone over V. The 
inner minimization is a nonconvex quartic optimization prob¬ 
lem in (zy) and is NP-hard in general. We relax the inner 
wtimization problem using the moment-relaxation approach 
jsj. Let a; = (l z'^ Re(L)^ Im(R)^). Let PolyVa;) 
denote the vector of all the monomials of degree upto i in (x) 
(with the first entry equal to the 0-degree monomial 1). For 
example Poly^ (x) = (l xi X2 ... x^ xiX2... ). 

Let m be the size of this Poly^ {x). We define a mo¬ 
ment vector y of the same size as Poly^ (x) and the linear 
operator on the space of all degree-4 polynomials in (x): 
Cy {EZlCXo\yt{x))=EZlCiyi. 


b -g 

g b 


the scaled Jacobian be- 
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We also define the localizing matrices (for i = 1,2): 
Xa = Poly' {x) (Poly' (a:))^. 

nmxmin tr (WCy (J^ {V) 

Subject to 

Cy ((Re(PiP/) - Xn) h 0 

Cy{{Re{ViV/)-^ij\Vj\^)X^i)yO 
((ILP -vf) Xii) = 0,i G pv 
Cy{{\Vo-lf)Xu)=0 
Cy {{z'^z - l) Xu) = 0 
Cy {X 21 ) h 0 

2/1 = 1, IIW^Lue < 1 

This is a convex-concave saddle point problem with compact 
feasible sets. Thus, we have strong duality, we can switch the 
min and max and reduce the problem to: 


min t 

y,t 

(14a) 

Subject to 

(14b) 

Sy(£j, {JF{V)zz^))^tI 

(14c) 

Cy ((Re iViVj*) - TiilLP) Xu) h 0 

(14d) 

Cy{{Re{ViVjn-lij\Vj\^)Xu)hO 

(14e) 

Cy{{\Vi\^-vyXu)=0,iGpv 

(14f) 

Cy{{\Vo-lf)Xu )=0 

(14g) 

Cy {{z'^z - l) Xu) = 0 

(14h) 

Cy {X21) h 0 

(14i) 

2/1 = 1 

(14j) 


If the above problem has a positive optimal value, then we 
obtain a certificate for ( |14c| ) is a semidefinite program in y 
that can be solved efficiently using interior point methods 0- 
W is given by the dual variable corresponding to the constraint 

( [Hel l. 

In section |V-B| we discuss ideas on scaling the approach to 
large networks by exploiting sparsity, and list the sizes of the 
resulting semidefinite programming problems. 


IV. Discussion 


The AC power fiow equations are fundamental to all aspects 
of power systems operations and planning, and several decades 
of research have gone into developing algorithms for solving 
the power flow problem. A discussion of the Newton’s method 
(by far the most popular algorithm developed for solving 
the PF equations) and its variants can be found in standard 
textbooks on power engineering p0| . The idea here is to 
update the power fiow variables according to: 

(z + 1) ^ (i) - 7Jt{JF { z)))“'f (z))‘) 


For a small enough step-size ry = if inverse Jacobian 
remains bounded ( (^Jp ^ < k for some /i: > 0), 

then the algorithm will converge to a power fiow solution. 
Given a good initial guess for Newton’s method con¬ 

verges rapidly. However, if the Jacobian becomes close to 


singular, then Newton’s method can behave badly. In general, 
Newton’s method can exhibit very complicated behavior and 
have fractal basins of attraction (TT) Several approaches 
(damping, trust region methods etc.) have been proposed and 
studied to improve the stability and convergence of Newton’s 
method. 

In the context of power systems, optimal multiplier methods 
were proposed |T^ GD to prevent divergence of Newton’s 
method. These adapt the choice of step length in Newton’s 
method to ensure that the algorithm does not diverge, although 
it may not converge to a power fiow solution even if the 
solution exists. A recent survey of the optimal multiplier 
algorithms along with numerical comparisons can be found 



An alternative approach to solving the PF equations is based 
on homotopy, where the power fiow problem is first solved 
for an “easy” injection vector and the injections are 

changed gradually while tracking the power fiow solution until 
the actual injections p^q are reached. This approach is called 
continuation power flow in the context of power systems 0- 
It was initially claimed that this algorithm is capable of finding 
all power fiow solutions, but this claim was shown to be 
incorrect recently fT^. N umerically this approach has been 
shown to be effective |17|, although theoretical guarantees are 
still lacking , to the best of our knowledge. 

Recently, a new approach based on complex analytic contin¬ 
uation techniques has generated a lot of interest and been 
shown to be effective for practical power systems problems. 
However, again, theoretical analysis of the conditions under 
which it is guaranteed to work are yet to be established.. 

A rigorous approach combining homotopy and tools from 
algebraic geometry was proposed 0 to find all power fiow 
solutions. However, computational scalability of this approach 
is currently limited to small networks. 

In l [2Q| , the authors propose a semidefinite programming 
(SDP) relaxation approach to solve the power fiow equations. 
The authors characterize the set of voltage vectors such that 
the SDP has a rank-1 solution and hence recovers a physically 
valid power fiow solution. This result is similar in spirit to 
ours. However, the set of recoverable voltages is specified by 
a nonconvex nonlinear matrix inequality which does not have 
a simple interpretation in terms of operational constraints on 
voltages. An exact comparison of the relative power of our 
approach and the approach proposed in this paper is a topic 
of future work. 

The closest works to what is described in the paper are 
1 [2T| and j^. In l [2T| , we studied the case of losses power 
networks Re(F^j) = 0 for every (i,jf) G £. In l [^ , we 
have shown that in the case of lossless networks, the power 
fiow equations can be solved by solving a convex optimization 
problem: Minimization of the energy function. The results of 
0 are a special case of the results in this paper. For lossless 
networks, the power flow operator F is the gradient of the 
energy function for the structure preserving model of power 
systems dynamics l [22| and Jp is the Hessian. In this case, the 
domain of monotonicity of the power fiow operator coincides 
with the domain of convexity of the energy function. We 
obtained a sufficient condition, expressed as a nonlinear but 
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convex matrix inequality in characterizing the set over 

which the energy function is convex (or equivalently, the power 
flow operator is monotone). However, for lossy networks, the 
monotonicity domain can no longer be characterized analyti¬ 
cally, which led us to the computational procedure described 
in section inrDi 

In Q, we described an alternate approach to computing 
monotonicity domains. The main difference between the work 
here and that presented in Q is the choice of coordinates 
used to compute the power flow Jacobian. In this paper, we 
worked with the log-polar coordinates Vi = exp -f-j 6 >i). 
The alternative , explored in Q, is to work with cartesian 
coordinates Vi = Vf + jV^. The choice of coordinates has 
been a subject of several numerical studies. In our experience, 
using log-polar coordinates has two main advantages: 

1) The size of the Jacobian is smaller, since the voltage 
magnitudes at PV buses are flxed and do not need to be solved 
for. This signiflcantly speeds up the computational procedure 
described in III-CI 

2) Our experiments show that the size of the monotonicity 
domain is signiflcantly larger. The reason for this phenomenon 
is not clear yet, and will be a topic of future investigation. 


V. Numerical Results 

In this section, we numerically validate our results. Section 


V-A describes the monotonicity domains computed for various 


test networks. Section discusses the sizes of the semidef- 
inite programming problems produced while solving ( p^ and 
section [\^ compares our approach to Newton’s method and 
the default solver from matpower | [23| , which is a variant of 
Newton’s method with some step-size control. 


A. Extent of Monotonicity Domains 

We first consider the 3-bus network plotted in figure 
All three buses are PV buses with voltage magnitudes set at 1 
p.u. In figure ( p^ , the blue region is the closed region whose 
boundary is defined by det {Jf {V)) = 0. Any convex domain 
of strict monotonicity domain cannot intersect the boundary 
of this region, since at the boundary the Jacobian is singular 
which means that Sy {WJf {V)) 0 for any W. Thus, any 

convex domain of strict monotonicity that contains the point 
V = 1 (all voltages equal to 1 p.u with 0 phase, the zero-load 
PF solution) must be contained inside the blue region. Our 
approach computes a bound of 7 = .08, which defines the red 
region. As one can see from the figure, this covers almost the 
entire blue region, so that our estimate of the monotonicity 
domain is nearly tight. 

For larger networks, we And the smallest uniform value 
7 > 0 such that the optimal value of is positive with 
yij = 7. The values are listed in table The first row lists 
the case number (matpower case file), the second row the 
minimum value of 7 for the network (as a uniform bound on all 
transmission lines) and the third row turns 7 into a bound on 
phase differences, assuming that the voltage magnitudes can 
vary between .9 and 1.1 p.u at all PQ buses. The number in 
the third row is a simplification of the actual constraint (which 
couples voltage magnitudes and phases). More generally, there 



(d) 3 bus monotonicity domain. Blue Region: 
True Monotonicity Domain. Red Region: Esti¬ 
mated Monotonicity Domain 



Fig. 1: 39 bus network: Tradeoff between bound on voltage 
ratios and phase differences 


is a tradeoff between the bound on voltage magnitudes ratios 
and phase differences. We plot this tradeoff for the 39 bus 
network in figure ( [2c| ). This shows that as voltage magnitudes 
are allowed to fluctuate more, phase differences need to 
kept within smaller limits to remain within the monotonicity 
domain. 


Case Number 

9 

14 

30 

39 

Min 7 

.31 

.34 

.41 

.52 

Max \6i — 6j 

67.7 deg 

65.4 deg 

59.9 deg 

50.54 deg 


TABLE I: Minimum 7 for different matpower test cases 
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Case number 

9 

14 

30 

39 

57 

Size of SDP 

78 

78 

105 

136 

190 


TABLE II: Size of PSD Constraints in sparse moment relax¬ 
ation 


B. Computation Time 

Solving the moment relaxation •El is the most expensive 
part of this approach. However, the situation is made better by 
the following factors: 

1) : This computation only depends on the network parameters 
and not on the specific injection profiles. Thus, the computa¬ 
tion can be done offline as long as the network topology is 
fixed. Further, even if the network topology changes (because 
of transmission switching, loss of a generator/line etc.), one 
can do offline computations for a large set of topology 
scenarios. 

2) The problem has a lot of sparsity. The Jacobian inherits 
the sparsity of the network, and each constraint only involves 
a small number of variables. One can exploit the sparsity 
in the moment relaxation to get more tractable Semidefinite 
Programs. Without sparsity, the size of the SDP would grow 
as where n is the number of buses in the network, 
which becomes intractable for n beyond 6 or so. However, 
by exploiting sparsity using the ideas described in p4| , we 
can reduce the size of the resulting SDPs significantly. We list 
the size of the largest PSD constraint for all the IEEE cases we 
tested in table One can see that the size grows gracefully. 

With the above considerations, we can solve the 39 bus 
network certification problem in about 20 minutes on a Mac- 
BookPro 2.6 GHz Intel Core i7 notebook. We believe that 
by using the following ideas, the approach can be scaled to 
networks with several thousand buses: 

1) In recent work pSj | |^ , the authors have shown that 
moment relaxations of the OPE problem can be solved for 
the networks with several thousand buses by applying higher 
order moment relaxations selectively. 

2) It is possible to consider more restrictive conditions than 
positive definiteness of Sy {WJp {V)). For example, one can 
require that this matrix is diagonally dominant or generalized 
diagonally dominant | [27| . This would give potentially more 
conservative results, but allow the SDP constraint to be re¬ 
placed by an LP (linear program) or SOCP (second order cone 
program), thus giving a more scalable approach. 

3) In the literature on graphical models and belief propagation 
| [28| has studied LP relaxations of sparse integer programs. By 
discretizing the space of voltages, similar techniques may be 
applied the problems studied here as well. 

C. Comparison to Newton's Method 

We also compare the performance of our method to a 
naive implementation of Newton-Raphson and the default 
matpower power fiow solver For a given network, we 
generate a random voltage phasor vector with voltages at the 
PQ buses sampled uniformly from the interval (.9,1.1) and 
phases at the non-slack buses sampled uniformly from the 
interval {—0,0) where 0 is a. bound on the absolute value 


of the voltage phase at each non-slack bus. We then form the 
injection vector corresponding to the voltage phasor and solve 
the resulting power fiow problem using 3 methods: a) The 
monotone operator method described in this paper, b) A naive 
implementation of Newton-Raphson with a constant unit step- 
size, initialized with all voltage phases equal to 0 and voltage 
magnitudes at PQ buses equal to 1 p.u, c) Matpower’s default 
PE solver. We generate a 100 random instances for each value 
of 0 and check whether each solver succeeded in finding a 
power fiow solution. We then plot the probability of success 
of each solver as a function of 0. The results for the 9-bus, 
IEEE 14-bus and IEEE 39-bus network (taking data from the 
matpower case files) are plotted in figure The results show 
that the monotone operator method is superior to the others. 
Further, it has the advantage of certifying the non-existence 
of a solution in the operationally relevant domain of voltage 
phasors, when it fails, a guarantee the other methods cannot 
provide. 

VI. Conclusions 

We have developed a new approach to solving the power 
fiow equations. The main feature of the new approach is its 
ability to find efficiently a solution within the “monotonicity 
domain” or prove that no solutions exist inside this domain. 
Our numerical results show that the monotonicity domain 
typically contains all power fiow solutions satisfying standard 
operational constraints on voltage magnitudes and phases. 
The computation of the monotonicity domain can be done 
offline and exploits the moment relaxation approach, which 
has recently been successfully applied to large scale optimal 
power fiow problems | [26| | [25| . 

Future work will focus on computational scalability of 
•El and extensions of this work to further applications: 
state estimation, optimal power fiow and small-signal stability 
analysis. In this regard, a related paper | [29| has recently been 
submitted, focusing on characterizing subsets of the feasible 
injection space in an OPE problem. 
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